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Abstract 

We study turbulence and Bose-Einstein condensation (BEC) within the two-dimensional Gross-Pitaevski 
(GP) model. In the present work, we compute decaying GP turbulence in order to establish whether BEC 
can occur without forcing and if there is an intensity threshold for this process. We use the wavenumber- 
frequency plots which allow us to clearly separate the condensate and the wave components and, therefore, to 
conclude if BEC is present. We observe that BEC in such a system happens even for very weakly nonlinear 
initial conditions without any visible threshold. BEC arises via a growing phase coherence due to anihilation 
of phase defects/vortices. We study this process by tracking of propagating vortex pairs. The pairs loose 
momentum by scattering the background sound, which results in gradual decrease of the distance between 
the vortices. Occasionally, vortex pairs collide with a third vortex thereby emitting sound, which can lead 
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to more sudden shrinking of the pairs. After the vortex anihilation the pulse propagates further as a dark 
sohton, and it eventually bursts creating a shock. 

1 Background and motivation 

For dilute gases with large energy occupation numbers the Bose- Einstein condensation (BEC) Pl 0] |S] can 
be described by the Gross-Pitaevsky (GP) equation m^]: 

i^-t + A^f- = 0, (1) 

where ^ is the condensate "wave function" . GP equation also describes light behavior in media with Kerr 
nonlinearities. 

Many interesting features were found in GP turbulence in both the nonlinear optics and BEC contexts 
ini El IHl El El E] ■ Initial fields, if weak, behave as wave turbulence (WT) where the main nonlinear process 
is a four-wave resonant interaction described by a four-wave kinetic equation 7 . This closure was used 
in [51IS1IH1 to describe the initial stage of BEC. It was also theoretically predicted that the four- wave WT 
closure will eventually fail due to emergence of a coherent condensate state which is uniform in space 9 . 
At this stage the nonlinear dynamics can be represented as interactions of small perturbations about the 
condensate state. Once again, one can use WT to describe such a system, but now the leading process 
will be a three- wave interaction of acoustic- like waves on the condensate background 9[. Coupling of such 
acoustic turbulence to the condensate was considered in | 15| which allowed to derived the asymptotic law of 
the condensate growth. 

In the stage of transition from the four- wave to the three- wave WT regimes, which itself is a strongly 
nonlinear process involving a gas of strongly nonlinear vortices, was studied. These vortices anihilate and 
their number reduces to zero in a finite time, marking a finite-time growth of the correlation length of the 
phase of ^' to infinity. This is similar to the Kibble-Zurek mechanism of the early Universe phase transitions 
which has been introduced originally in cosmology ^1E|- It has been established that the vortex anihilation 
process is aided by the presence of sound and it becomes incomplete if sound is dissipated. Fourier transforms 
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in both space and time were analysed using the wavenumber-frequency plots which, in case of weak wave 
turbulence, are narrowly concentrated around the linear dispersion relation lo ~ ujk- At the initial stage, 
narrow (fc, ct;)-distributions around w = cjfe = fc^, were seen whereas at late evolution stages we saw two narrow 
components: a condensate at horizontal line lu = (p) = and an acoustic component in proximity of 

the Bogolyubov curve Qk = (p) + v^fc^ + 2{p)k'^. 

In Pni, the system was continuously forced at either large or small scales because this is a classical WT 
setting. WT predictions were confirmed for the energy spectra of GP turbulence. However, it remained 
unclear if presence of forcing is essential for complete BEG process, and whether there is any intensity 
threshold for this process in the absence of forcing. This questions are nontrivial because, in principle, even a 
weakly forced system could behave very differently from the forced one due to an infinite supply of particles 
over long time. 

In the present paper we will examine these questions via numerical simulations of the 2D GP model 
without forcing. In addition, we will carefully examine and describe the essential stages of the typical route 
to the vortex anihilation leading to BEG. In many ways our work is closely related with numerical studies of 
decaying 3D GP turbulence of Ref . fTSI , where transition from the 4- wave WT regime to condensation and 
superfluidity was studied by visualising emergence and decay of a superfluid vortex tangle. However, our 
work addresses additional issues which were not discussed in the previous papers, particularly emergence and 
dispersive properties of the 3-wave acoustic turbulence, crucial role of sound for the vortex decay process, 
separating the wave and the condensate components using novel numerical diagnostics based on {k,Lu)- 
distributions. 

Below, we will only describe our numerical setup and results. For a summary of WT theory and its 
predictions in the GP context we reffer to our previous paper jlUj . 

2 Setup for numerical experiments 

In this paper we consider a setup corresponding to homogeneous turbulence and, therefore, we ignore finite- 
size effects due to magnetic trapping in BEG or to the finite beam radii in optical experiments. For numerical 
simulations, we have used a standard pseudo-spectral method 18 for the 2D equation (Q: the nonlinear 
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term is computed in physical space while the linear part is solved exactly in Fourier space. The integration 
in time is performed using a second-order Runge-Kutta method. The number of grid points in physical space 
was set to N X N with N — 256. Resolution in Fourier space was Ak — I-k jN . Sink at high wave numbers 
was provided by adding to the right hand side of equation 1^ the hyper-viscosity term V^)"?/'- Values of 
V and n were selected in order to localized as much as possible dissipation to high wave numbers but avoiding 
at the same time the bottleneck effect, - a numerical artifact of spectrum pileup at the smallest scales 
Note that importance of introducing the small-scale dissipation to eliminate the bottleneck effect has long 
been realised in numerical simulations of classical Navier-Stokes fluids, and it was also recently realised in 
the context of GP turbulence in 14 . We have found, after a number of trials, that v — 10^^ and n — % 
were good choices for our purposes. Time step for integration was depended on the initial conditions. For 
strong nonlinearity smaller time step were required. Numerical simulations were performed on a PowerPC 
G5, 2.7 Ghz. Initial conditions were provided by the following: 

V'(fc.,fc„t-0) = ^4=e-^e'* (2) 

V TT-^/^lT 

where fc = y^/c^ ~^^'y- is a real number which was varied in order to change the nonlinearity of the initial 
condition for different simulations. (\) ~ 4>{kx, ky) are uniformly distributed random numbers in the interval 
[0, 2tt]. The simulations that will be presented here have been obtained with kQ~ 35 Afc and a = 5Afc. The 
nonlinearity of the initial condition was measured as e=fcoa. We have performed simulations ranging from 
Emm = 0.018 to Emax = 1-3, SO wc havc Spanned almost 2 decades of nonlinear parameter in the initial 
conditions. 

3 Numerical results 
3.1 Evolving spectra 

We start by examining the most popular turbulence object, the spectrum. 
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Figure 1: Spectra at different times for the initial condition characterized by e = 1.3 

Since we study the setup corresponding to BEC, we start with a spectrum concentrated at high wavenumbers 
leaving a range of smaller wavenumbers initially empty so that it could be filled during the evolution. In 
figures n and |5] we show the spectra at different times for the strongest and the weakest nonlinear initial 
condition we have analyzed. 

At early stages, for both small and large initial intensities, we see propagation of the spectrum toward 
lower wavenumbers. However, we do not observe formation of a scaling range corresponding to the wave 
energy equipartition rik ^ 1/ujk — l/k'^ (i.e. energy density tOkUk is constant in the 2D /c-space) as it was the 
case in the simulations with continuous forcing |10| . 

At later stages, no matter how small the initial intensity is, the low-k front reaches the smallest wavenum- 
ber and we observe steepening reaching slope '-^ — 3.5 for large initial intensities (figure^ and ^ — 1 for the 
weakest initial data (figure [21). This corresponds to WT breakdown and onset of BEC. However, the infor- 
mation contained in spectrum Uk is very incomplete as it does not allow to distinguish between the coherent 
condensate and random waves that may both occupy the same wavenumber range. Thus, we turn to study 
the direct measures of condensation such as the correlation length and the wavenumber-frequency plots. 
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Figure 2: Spectra at different times for the initial condition characterized by e = 0.018 
3.2 Explosive growth of correlation length 

By definition, condensate is a coherent structure whose correlation length is of the same order as the bounding 
box. We define the correlation length directly based on the auto-correlation function of field 

C^(r) = (7^*(x)7^*(x + r))/(7^^'(x)2), (3) 

where TZ'^ denotes the real part of ^' (result based on the imaginary part would be equivalent). Correlation 
length A can be defined as 

pro 

/ C^{r)dr, (4) 
Jo 

where tq is the first zero of C^(r). Note that initially C^{r) can strongly oscillate, which is a signature of 
weakly nonlinear waves. However, only one oscillation (i.e. within the first zero crossing) is relevant to the 
condensate, which explains our definition of A. Figure O shows evolution of 1/A^ which, as we see that A 
always reaches the box size wich is a signature of BEC. 
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Figure 3: Evolution in time of the correlation length for different nonlinearity. 
3.3 Wavenumber-Frequency plots 

As we discussed above, the spectra cannot distinguish between random waves and coherent structures es- 
pecially when they are present simultaneously and overlap in the fc-space. Besides, the spectra do not tell 
us if the wave component is weakly or strongly nonlinear. To resolve these ambiguities, following |10| . let 
us perform an additional Fourier transform over a window of time and examine the resulting (fc, a-')-plots of 
space-time Fourier coefficients. 

Figure 21 corresponds to an early time of the system with relatively weak initial intensity (e = 0.4). We 
see that the distribution is narrowly concentrated near lo — k"^ which indicates that these waves are weakly 
nonlinear. The weak nonlinear effects manifestate themselves in a small up-shift and broadening of the (fc, lo)- 
distribution with respect to the lo = kF' curve. For sufficiently small initial intensities, these early stages of 
evolution are characterised by weak 4-wave turbulence. The breakdown of the uj = k'^ curve at high fc's 
occurs due to the the numerical dissipation in the region close to the maximal wavenumber (this component 
is weak but clearly visible because the color map is normalised to the maximal value of the spectrum at each 
fixed k). 
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Figure 4: (A;,cj) plot for the initial stage for the case of e = 0.4. Solid curve shows the Bogolyubov dispersion 
relation. 




Figure 5: (/c, uj) plot for the late stage for the case of e = 0.4. Solid curve shows the Bogolyubov dispersion 
relation. 




Figure 6: (w — k) plot for the latest stages for the most nonlinear case, e = 1.3. Solid curve shows the 
Bogolyubov dispersion relation. 

Figure shows a late-time plot for the same run (i.e. e — 0.4). The late-stage (fc, Ci;)-plots for the most 
nonlinear intensities are in shown in figure El We see that in both cases we now see two clearly separated 
components quite narrowly concentrated around the following curves: 

• (A) A horizontal line with w « (p) , 

• (B) The upper curve which follows the Bogolyubov curve w = fifc = (p) + -y/fc^ + 2{p)k?. 

Component (A) corresponds to BEC. Its coherency can be seen in the fact that the frequency of different 
wavenumbers is the same. Note that usually BEC is depicted as a component with the lowest possible 
wavenumber in the system, whereas in our case we see a spread over, although small, but finite range of 
wavenumbers. This wavenumber spread is caused by few remaining deffects/ vortices. 

Curve (B) corresponds the Bogolyubov sound-like waves. We see that the the (fc, a;)-distribution is quite 
narrow and close to the Bogolyubov dispersion curve, which indicates that these waves are weakly nonlinear. 
However, now these weakly nonlinear waves travel on a strongly nonlinear BEC background. This is a 
three- wave acoustic weak turbulence regime P 1151 ITT)] . 
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Figure 7: Total wave spectrum (solid line), condensate wave spectrum (dotted line) and background wave 
spectrum (dotted line) 

3.4 Separating the condensate and the wave components 

Using the Wavenumber- Frequency plots we can separate BEC and the wave component and plot their spectra 
separately. In figure [7| we show the spectrum for the case of e = 1.3 at late time for the condensate and the 
rest of the wave field. These spectra have been obtained by integrating the {u; — k) from to a threshold cjc 
(this corresponds to condensate) and from oj^ to the maximum value of oj considered. In the present case, in 
order to separate the condensate, lUc was set to 1.5 (see figure 15)). As is clear from the the figure, most of the 
energy is concentrated in the condensate. 

3.5 Typical events leading to vortex anihilation 

Vortices ar phase deffects of and, therefore, the correlation length growth is intimately connected with the 
decrease of the total number of vortices. As argued in 10 that Nyortices ^ To give an illustration of 

the anihilation process we show in figure|Sltwo snapshots of \'i>{x,y,t)\ at different times. The vortices are 
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Figure 8: \^{x,y,t))\ at two different instants of time 



seen as blue spots in these snapshots and we see a considerable decrease of their number at the later time 
(on the right). It was shown in 10 that Bogolyubov sound is an essential mediator in the vortex anihilation 
process, and that introducing a sound absorption can lead to frustration and incompleteness of the BEC 
process. 

Let us now examine in detail the typical sequence of events leading to the vortex annihilation by tracking 
vortex pairs that are destined to anihilate. The vortex pair motion is best seen in a computer generated 
movie which is available upon request from the authors. In Figures 1^ IIUI and 1111 wc show a representative 
sequence of frames from this movie. 

We can see that the vortex pair forms in frames a to c so that the distance between the vortices in the 
pair is considerably less than distance to the other vortices. Such a pair propagates like a vortex dipole in 
fluid, as seen in frames c to f. During this propagation, the vortex pair scatters the ambient sound waves 
thereby transferring its momentum to the acoustic field. This momentum loss make vortices get closer to 
each other and, therefore, move faster as a pair. This process can be interpreted as a "friction" between the 
vortices and a "normal component" (phonons). It is easy to show that such process leads to the change of 
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Figure 13: Distance between a pair vortices as a function of time for the vortex pair and the time range of 
Figure M 



distance d(t) between the vortices hke 



d(t) = a\/t^ - t, 



(5) 



where coefficient a is proportional to the energy density of sound and t"^ is the anihilation time. Figure [T^ 
shows evolution of the inter- vortex distance d{t) , calculated for the vortex pair and the time range of Figure 
El We see that for the time range when the vortex pair is more or less isolated from the other vortices, 
11 < i < 40, the inter- vortex distance shrinks in qualitative agreement with law (jSJ. However, immediately 
after that, in frame h in Figure IIUI the vortex pair collide with a third vortex and suddenly shrink and 
anihilate. Further on the movie we see that the vortex pair momentum was not completely lost and it keep 
propagating a Jones- Roberts dark soliton |19| . In between of frames h and o the soliton "cannot make up 
his mind" oscillating between the state with and without vortices near the vortex anihilation threshold, until 
frame p where it re-emerges as a vortex pair. It collides with yet another vortex and changes its propagation 
direction to 90 degree in between of frames q and r, it annihilates again in frame s, propagates as a dark 
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soliton until frame v. Eventually, this soliton is weakened due to further sound generation and scattering 
and it becomes too weak to maintain its stability and integrity. At this point it bursts thereby generating a 
shock wave as seen in frames v to y. 

Summarizing, we can identify the following important events on the route to vortex anihilation: 

• Gradual shrinking of the inter-vortex distance due to the sound scattering when the vortex pair is 
sufficiently isolated from the rest of the vortices, 

• Sudden shrinking events due to collisions with a third vortex and resulting sound generation (similar 
three- vortex event was described in SO;), 

• Post-anihilation propagation of dark solitons, 

• Occasional recovery of vortex pairs in dark solitons which are close to the critical amplitude, 

• Weakening and loss of stability of the dark solitons resulting in a shock wave. 

4 Conclusions 

We have computed decaying GP turbulence in a 2D periodic box with initial spectrum occupying the small- 
scale range. We observed that BEG at large scales arises for any initial intensity without a visible threshold, 
even for very weakly nonlinear initial conditions. BEG was detected by analysing the spectra and the 
correlation length and, most clearly, by analysing the {k, us) plots. On these plots BEG and Bogolyubov waves 
are seen as two clearly distinct components: (i) BEG coherently oscillating at nearly constant frequency and 
(ii) weakly nonlinear waves closely following the Bogolyubov dispersion. By separating BEG and the waves 
we observed that most of the energy at late times is residing in the BEG component at most of the important 
scales except for the smallest ones. 

We also analysed the typical events on the path to vortex anihilation, - an essential mechanism of BEG. 
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